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Abstract 

We study a multi-determinant approach to tiie time evolution of the nu- 
clear wave functions (TDMD). We employ the Dirac variational principle 
and use as anzatz for the nuclear wave-function a linear combination of 
Slater determinants and derive the equations of motion. We demonstrate 
explicitly that the norm of the wave function and the energy are conserved 
during the time evolution. This approach is a direct generalization of the 
time dependent Hartree-Fock method. We apply this approach to a case 
study of ^Li using the N3L0 interaction renormalized to 4 major harmonic 
oscillator shells. We solve the TDMD equations of motion using Krylov sub- 
space methods of Lanczos type. We discuss as an application the isoscalar 
monopole strength function. 
Pacs numbers: 21.60.-n, 24.10.Cn, 31.70.Hq 
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1 Introduction. 



The time dependent Hartree-Fock method (TDHF) and its quasi-particle gener- 
alization, the time dependent Hartree-Fock-Bogoliubov method (TDHFB), are 
central tools in studying nuclear dynamics (see for example ref. [1]), whereby 
the nuclear wave function is approximated by a single Slater determinant or by 
a quasi-particle determinant wave-function. In the small amplitude limit they re- 
duce to the random phase approximation (RPA) and to the quasi-particle random 
phase approximation (QRPA) which are commonly used to describe nuclear exci- 
tation, as for example giant resonances. 

Despite the enormous matrix dimensions, these limits are nowadays solved us- 
ing efficient Krylov projection techniques of Arnoldi type (see for example ref. [2] 
for recent applications). Recently the time dependent coupled-cluster method 
(refs. [3], [4]) has been revisited and applied to light nuclei using a similarity 
group transformed N3L0 interaction (ref. [5]). 

In this work we discuss a Time Dependent Multi-Determinant (TDMD) ap- 
proach whereby the nuclear wave function is approximated by a linear combina- 
tion of Slater determinants. This approach is the time dependent version of the 
Hybrid Multi-Determinant (HMD) approach (refs. [6], [7]). To the author knowl- 
edge, this approach has never been considered in nuclear physics. In this sense 
this is an exploratory study. Our starting point is the Dirac variational principle 
which, as well known, leads to the time-dependent Schrodinger equation or the 
time-dependent equations of motion as the TDHF (ref. [8]) within the approxima- 
tion of the nuclear wave function with a single Slater determinant. Within the 
Dirac variational principle we derive the equations of motion and demonstrate ex- 
plicitly that the time evolution conserves the norm of the wave function and the 
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energy. Since the equations of motion are, as we shall see, of the type iLtp — R 
where R is an energy gradient ^ is a sum of Slater determinants and L is a ma- 
trix (of rather large dimensions) embodying the time-derivative of the norm of the 
wave-function (which will be discussed in detail below), the actual evaluation of 
the wave function as a function of time is performed using the Direct Lanczos 
method (DL) for the solution of a rather large linear system. The DL method is 
a Krylov subspace type of method to solve linear systems (an excellent review of 
these methods can be found for example in ref. [9]). The eigenvalue version of 
these methods are the familiar Lanczos method used in the shell model approach 
to nuclear structure and the Arnoldi method used in solving the RPA or QRPA 
eigenvalue problems. The basic idea of these methods is that although we cannot 
store a matrix (e.g. the nuclear Hamiltonian matrix) we can evaluate the matrix to 
vector product easily. Although L is not as large as the shell model Hamiltonian 
matrix it can hardly be stored except in very simple cases. However the matrix 
to vector product appearing in the equations of motion is trivial to evaluate, and 
therefore the Laczos method is the ideal one. We solve the equations of motion, as 
an exploratory study, in the case of ^Li using the N3LO interaction renormalized 
to 4 major oscillator shells with the Lee-Suzuki (ref.[10],[l 1]) method, in order to 
reduce the otherwise very large single-particle space. We use the time-dependent 
wave function thus obtained to evaluate strength functions. Our ultimate goal is 
to extend ab-initio methods to time-dependent problems, such as the evaluation of 
strength functions, starting from the nucleon-nucleon interaction. 

The outline of this paper is as follows. In section 2 we derive the equations 
of motion using the Dirac variational principle and prove that these equations of 
motion conserve the norm and the energy of the nuclear wave function. In section 
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3 we discuss the numerical method and in section 4 we discuss the appUcation to 
the nuclear strength function using the boost method. 



2 The time dependent variational principle. 

The Dirac time-dependent variational principle states that the time evolution of 
the nuclear wave function is obtained by varying the action 

Si = dtCi = /*' dt[ih <ip\tp> -< ip\H\ip >] (la) 
Jtl Jtl 

or equivalently 

,^2 = /*' dt£2 = r dt[-ih <ij\ij> - < ij\H\ij >] (16) 
Jtl Jtl 

with respect to |^ > and < ^\ independently, under the constraint that the wave 
function is held fixed at the initial and final times ti and ^2- The most general 
wave function in the Hilbert space will then give the Schrodinger equation and its 
complex conjugate. The TDHF approximation is obtained if the wave function 
is assumed to be a Slater determinant. In what follows we drop h with the un- 
derstanding that the unit of time is lMeV~^ ~ 6.6 x 10~^^sec. We consider the 
ansatz 

\ij >= ^ \Us > (2) 

5=1 

where \ Us > is a Slater determinant and A^^ is the total number (space dimension). 

These Slater determinants are of the most generic type and are written as 

\Us >= cIscIs..Jas\0 > (3) 
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A being the number of particles, S labels the Slater determinant and 

cL= E Ui,asal (a = 1,2,.., A) (4) 

in the above equation, a| is the creation operator in the single-particle (harmonic 
oscillator) state i, Ngis the number of the single-particle states and U is the single- 
particle wave-function in the h.o. representation. We assume that neutrons and 
protons are not mixed, i.e that each Slater determinant is a product of a neutron and 
a proton Slater determinant. The ansatz (3)-(4) is the same as in the Hybrid multi- 
determinant method (refs.[6],[7]). Here we omit the projector to good quantum 
numbers (angular momentum and parity) solely for clarity. In principle we could 
apply the variational principle including projectors to good quantum numbers. 
The Slater determinants are not orthogonal to each other and are 'deformed', that 
is they do not have good quantum numbers. At the initial time they could be 
the result of a partially converged variational calculation as given by the HMD 
method. The coefficients of the linear combination in eq.(3) are set to 1 since they 
can be included in the definition of the Ui^aS- 

We seek to find the equations of motion for the single-particle wave functions 
U. In what follows, the complex conjugate of the Slater determinant \Us > is 
denoted as < 0\cas'---Cis', where 

Ca,S' = ^Va,iS'ai (5) 
i 

and V is the hermitian conjugate of the matrix U. We do this in order to use simple 
matrix notations. Hence the indices of U are i, a, 5" and the indices of V are 
a,i,S'. In the following greek letters will label particles and latin letters single- 
particle states. The symbols S and S' will label different Slater determinants. The 
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Dirac variational principle gives (the bra will be denoted as < y | ) 

il2^viS') < Vs>ps >= Sv{S')J2< Vs\H\Us > (5a) 
-iJ2^u(s) < Vs'\Us >= Suis)J2< ^s\H\Us > (56) 

S' S' 

where we we shown explicitly the quantities which are varied. In what follows we 
quote the results for the overlaps and for the matrix elements of the Hamiltonian 
(of. ref.[6]). The Hamiltonian is 

H ^^Yl Hijkia\a]aiak (6) 

ijkl 

where we recast the one-body term into the two-body interaction, as done in shell 
model calculations and the matrix elements H are antisymmetrized (i.e. Hij^i = 
—Hijik). For any V and U, (relative to the Slater determinants S' and S) let us 
define 

G^(VU)-\ W^GV, X^UG, p^UGV, F^l-p (7) 

The matrix G has indices a,/3 — 1, 2, .., A. The matrix W has indices a, i, the 
matrix X has indices i , a while p and F have indices i,j = 1,2, ..,Ns. The matrix 
p is the generalization of the density matrix in TDHF and satisfies the relations 
trp = A and p^ — p for any S' and S, as it can easily be verified. Then (cf. 
(ref.[6])) we have 

< V\U >= det{VU) (8) 
< V\H\U >=< V\U > tr(rp) (9) 
where the matrix F is given by 

Fjj = ^ HpiqjPqp (10) 
pq 
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Let us note that the exchange tenn is the same of the direct since the matrix el- 
ements are antisymmetrized. Let us note that the arrays in eqs.(7)-(10) are eval- 
uated with Vs' and Us. Also let us note that the single particle wave functions 
are not orthogonal to each other. Only in the case of a single Slater determinant 
they can be made orthogonal. The equations of motion eqs.(5a),(5b) (EOMl and 
E0M2) can now be derived using the matrix identity 

Sdet{M) = det(M)tr(M-^5M) (11) 

Then it is easy to verify that 

< V\U >=< V\U > tr(GVU) (12a) 

< V\U >=< V\U > tr(GVU) (126) 
and that the expUcit form for EOMl is (using the identity SM~^ = -M~'^SMM~^) 

i ^ det{VU){Xi^W^r + FirG^^)Ur^,s = Edet(V^t/)(Xia^ + 2(FrX),, (13) 

riJ.,S S 

where S is the energy functional 

S = tr(rp) (14) 
The equation of motion E0M2 can be obtained in the same way and is given by 

-I ^ det{VU){Wa^Xr^ + G^|_,Fr^)V^rS' = Yld<^U){SWa^ + 2{WTF)^,) 
H,r,S' S 

(15) 

These equations need a few comments. First, if we recast then in a schematic 
matrix notation 

iL^'^il = (16a) 
-iL(2)T/ = (166) 



the dimension of the hnear systems to be solved can be rather large. For example 
In the case of '^'^Mg with 7 major shells (iY^ = 168) for 10 Slater determinants, 
the matrix L is 20160 x 20160 (for both neutrons and protons), for many more 
major shells or heavier nuclei, the storage of this array in the computer RAM is 
a problem. Moreover these matrices seem to have a separable structure. L^^^ for 
example contains a separable term in the (ia) [jir) indices and another term sep- 
arable in the {ir) {jia) indices. This implies that although we may not be able to 
store the matrix L, we can very easily perform the matrix to vector product. We 
only needs to store the matrices X, W, F and G, in the above Mg case of dimen- 
sion 168 X 12, 12 X 168 and 168 x 168 for every Slater determinant. In the past 
few decades, linear systems of this kind, where the matrix cannot be stored but we 
can easily compute matrix to vector product have received a major attention in ap- 
plied mathematics using the so called Krylov subspace techniques. It is precisely 
the same kind of techniques one uses in standard shell model calculations. They 
will be briefly recalled in the next section. A systematic treatment can be found 
in ref. [9] (note however that in ref [9] the convention for the scalar product is 
< x\y >—Y, XiUi). Equations of motions EOMl and EOM2 are equivalent. The 
matrix L^^^ and L^^^ are hermitian. 

One can show that the norm of the wave function is preserved during time 
evolution. From eqs. (8) and (11) one has 

d<'il^\ilj> /dt = Y,< y\U > tr[G(VU + VU)] (17) 

SS' 

with the understanding that S' refers to V and S to U. From EOMl eq. (13), 
multiplying by VaiS' and summing over the indices one has 

iJ2<V\U> [tr(p)tr(GVU) + tr(FUGV)] = 

S'S 

8 



Y,<y\U> Mp)S + 2tr(Frp)] (18) 

S'S 

From EOM2 of eq. (15), multiplying by Uias and summing over the indices one 
has 

-iY,<y\U> [tr(p)tr(GVU)+tr(XVF)] = ^ < V|U > [tr(p)£:+2tr(prF)] 

SS' S'S 

(19) 

Subtracting eqs. (17) and (19), since for any SS' tr(p) = A, 
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^ < V\U > [A tr[G(VU + VU)] + tr(FUW + XVF)] 



SS' 

Y,<V\U> 2tr(Frp - pFF) (20) 

S'S 

The right hand side of this equation is since F = 1 — p. Next, since tr(p) = 
for any S'S using the definition of p given in eq.(7) and the cyclic property of the 
trace, we have 

tr[G(VFU + VFU)] = (20) 

This equation can also be verified directly using the definitions in eq. (7). Hence 
from eq.(20), using the cyclic property of the trace, one has 

< y|C/ > tr[G(VU + VU)] = (21) 

SS' 

which is the time derivative of the norm (cf. eq.(17)). Next we shall prove that the 
energy is constant during the time evolution. We need to prove that 

d<'ijj\H\'ijj> /dt^O (22) 

since the norm of the wave function is a constant. Let us set 

n[v, U] =< VI^IV' >, o[v, U] =< VIV' > (23) 
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The Lagrangian associated to EOMl can be rewritten as 



A = iE^^a-^ (24) 



where a — (iaS) for brevity. EOMl can then be written as 

'\dv,du!^^- dv, ^^^^ 

for all h — {/3jS'). Similarly the Lagrangian associated with EOM2 can be recast 



as 



A—'EaT^H-^i (26) 



and EOM2 can be recast as 

-E^H^^ (27) 

Multiplying eq.(25) by Vf, and summing over the indices and multiplying eq.(27) 
by Ua and summing over indices and subtracting the two results we obtain 

which is precisely the time derivative of "H. These two constants of motion are 
a valuable test to check whether the equations of motion have been integrated to 
reasonable accuracy. 



3 A brief description of the numerical method. 



We solve numerically EOM2 (eq.l5) and eq.(16b)). The reason behind this choice 
is that the subroutines which are needed can be extracted from HMD computer 
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programs, which have been tested to very high accuracy. As pointed out in the 
previous section, it is not advisable to store the matrix L2. However we can easily 
evaluate L2V where v is any vector. Actually we can easily evaluate any power of 
L2 applied to v. The linear system of eq.(16b) can then be solved by projecting 
eq.(16b) into the subspace generated by the vectors v, L2V, (L2)^f , . . where v is an 
arbitrary trial solution of the linear system (known as Krylov subspace), followed 
by Graham-Schmidt orthonormalization. Since L2 is hermitian the projection of 
L2 gives a tridiagonal matrix (just as in the standard shell model method) and the 
linear system can then be efficiently solved. We have implemented the so called 
direct Lanczos method, the full detail of which (including the algorithm) can be 
found in ref. [9]. In our computer program the iterations stop when the residual 
vector —iL2V — R2 has a norm less than 10~^^. In this work we considered ^Li 
with the interaction given by the N3L0 NN potential renormalized using the Lee- 
Suzuki method to 4 major harmonic oscillator shells. Once V has been determined 
we solve the differential equation in time using a rank-4 Runge-Kutta method. The 
time step is typically IQ^^MeV^^ . In fig. 1 we show the errors in the energies 
and in the overlaps DE = E{t) - E{Q) and DO = ln(C>(i)/C>(0)) as a function 
of time. In fig. 1 we took snapshots every 500 (typically) time steps. Typically 
we ran the time evolution up to few hundreds MeV~^. The number of Lanczos 
iterations needed for convergence depends on the number of Slater determinants. 
For 1 Slater determinant (TDHF) we need about 4 Lanczos iterations to reach 
machine accuracy. This number increases as we increases the number of Salter 
determinants. For example for 3 Slater determinants we need typically 36 Lanczos 
iterations and for 5 Slater determinants we need about 50 iterations. In all runs we 
considered hio — 12MeV and we added to the Hamiltonian the standard center of 
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Figure 1 : Variation of the energy (in MeV) and variation of the norm as a function 
of time. We took 3 Slater determinants with 4 major shells at hcu = 12MeV 
Snapshots are taken every 0.5MeV~^. 
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mass Hamiltonian with a strength equal to 1. 



4 Strength functions. 

We evaluate the strength function for a one body operator Q 

S{E) = Y.\<n\Q\Q>\'5{E-El) 

n 

E* being the excitation energy of the n-th approximate eigenstate, with the boost 
method, as follows. First we determine the approximate ground-state of the sys- 
tem |0 >, then at time i = 0+ we boost the system with the unitary operator 

|V'(0+) =exp(i?7Q)|0> 

for sufficiently small rj. We then evolve this wave function in time by solving the 
equations of motion and evaluate 

Q{t) ^<m\Qm)> 

The strength function can be obtained using the Fourier transform of Q{t) (see for 
example ref. [12]) 

Jo 

for sufficiently large T such that e"'^^ is negligible via the relation 

S(E) = —Im(Q(E)) 
rjTT 

The width F is sufficiently small and such that very high frequency oscillations 
are smoothed out. We take typically F = O.lMeV. We need to evolve the system 

13 




Figure 2: Monopole strength function for Nyj = 1 (TDHF). 
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Figure 4: Monopole strength function for A^^ = 2 Slater determinants. 



16 



10 15 20 25 30 

E (MeV) 



Figure 5: Magnification of fig. 4. 
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Figure 6: Strength function for iVu, = 3 Slater determinants. 
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Figure 7: Magnification of fig. 6. 
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Figure 8: Strength function for Nyj = 5 Slater determinants. 

after the boost for about T = 100MeV~^. In this exploratory calculation we have 
used the isoscalar monopole operator Q — r'^. 

In figs. (2),(3) we show the results obtained in the case 1 Slater determinant 
(TDHF). Fig (3) is a magnification of fig (2). The strength function is dominated 
by the dominant peak at £^ ~ 17.8MeV. Some weaker peaks can be seen at 
E ~ O.lMeV, E ~ 7.56Mey and E ^ 27.6MeV. Using 2 Slater determinants 
we obtained the results shown in figs. (4) and (5).The dominant peak is now at 
E ~ 16.2MeV. The secondary maxima are atE ~ O.lMeV, 7.5MeV, 27.SMeV, 
almost on the same position of the TDHF case. Using 3 Slater determinants, we 
obtained the results of figs, (6) and (7). The main peak at £■ ~ 16.7 MeV shows 
considerable fragmentation around 15 MeV while the secondary peak at 7.5MeV 
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Figure 9: Magnification of fig. 8. 
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is nearly unchanged. The peak around 27 MeV has nearly disappeared and has 
moved to lower excitation energies. We considered in this work a maximum of 
5 Slater determinants, and the results for the strength function are shown in fig. 
8 and the details in fig. 9. The main peak is located at lY.lMeV closer to the 
TDHF case. However secondary peaks are at slight different energies and the 
strength shows a more pronounced fragmentation. 

We have not studied the strength function for a larger number of major shells 
and as a function of the h.o. frequency for increasing number of Slater determi- 
nants. Such a study is necessary in order to promote the method as an ab-initio 
method. Nor we have studied the small amplitude limit of the equations of motion 
in the case of many Slater determinants which is the equivalent of the RPA with 
several Slater determinants. Our purpose in this work is to define the method, 
solve the equations of motion and verify our computer programs. Such studies 
will be presented in future works. 
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